Loading the coordinates data

We need the sf package:

> if (! "sf" %in% rownames(installed.packages())) install.packages("sf")
> library(sf)
Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3

The coordinates of the houses of the Hanam cohort are in a zipping shapefile that can be downloaded as so:

> download.file("https://www.dropbox.com/s/yyg6gx0nycpv60b/HH%20point%20Hanam.zip?raw=1", "hanam.zip")

Once downloaded, unzip the file:

> unzip("hanam.zip")

Once unzipped read the shapefile:

> hanam <- sf::st_read("Household_point.shp")
Reading layer `Household_point' from data source `/Users/choisy/Dropbox/aaa/projects/Hanam/Household_point.shp' using driver `ESRI Shapefile'
Simple feature collection with 263 features and 22 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

And clean the disk

> file.remove(grep("^Household", dir(), value = TRUE))
[1] TRUE TRUE TRUE TRUE TRUE TRUE

and

> file.remove("hanam.zip")
[1] TRUE

Here the data, which is of sf class

> hanam
Simple feature collection with 263 features and 22 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
     Header Name          Descriptio          Type             Position     tempX        X        Y Altitude Depth Proximity Temperatur
1  Waypoint H001 14-JAN-09 9:20:51AM User Waypoint N20.49189 E105.92737 N20.49189 20.49189 105.9274    314 m  <NA>      <NA>       <NA>
2  Waypoint H002 14-JAN-09 9:18:10AM User Waypoint N20.49196 E105.92861 N20.49196 20.49196 105.9286    318 m  <NA>      <NA>       <NA>
3  Waypoint H003                <NA> User Waypoint N20.49663 E105.93138 N20.49663 20.49663 105.9314      0 m  <NA>      <NA>       <NA>
4  Waypoint H004 14-JAN-09 9:20:08AM User Waypoint N20.49191 E105.92738 N20.49191 20.49191 105.9274    317 m  <NA>      <NA>       <NA>
5  Waypoint H005                <NA> User Waypoint N20.49529 E105.93026 N20.49529 20.49529 105.9303      6 m  <NA>      <NA>       <NA>
6  Waypoint H006                <NA> User Waypoint N20.49530 E105.93292 N20.49530 20.49530 105.9329     -3 m  <NA>      <NA>       <NA>
7  Waypoint H007 14-JAN-09 2:10:03PM User Waypoint N20.49722 E105.92956 N20.49722 20.49722 105.9296      1 m  <NA>      <NA>       <NA>
8  Waypoint H009 14-JAN-09 3:03:29PM User Waypoint N20.49759 E105.92656 N20.49759 20.49759 105.9266      1 m  <NA>      <NA>       <NA>
9  Waypoint H010                <NA> User Waypoint N20.49672 E105.92964 N20.49672 20.49672 105.9296     -3 m  <NA>      <NA>       <NA>
10 Waypoint H011                <NA> User Waypoint N20.49684 E105.93064 N20.49684 20.49684 105.9306      2 m  <NA>      <NA>       <NA>
      Display_Mo Color    Symbol Facility City State Country Date_Modif Link Categories                  geometry
1  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9274 20.49189)
2  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9286 20.49196)
3  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9314 20.49663)
4  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9274 20.49191)
5  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9303 20.49529)
6  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA>  POINT (105.9329 20.4953)
7  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9296 20.49722)
8  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9266 20.49759)
9  Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9296 20.49672)
10 Symbol & Name Black Residence     <NA> <NA>  <NA>    <NA>       <NA> <NA>       <NA> POINT (105.9306 20.49684)

The coordinates of the houses can be extracted from that object and converted into a data frame as so:

> houses <- as.data.frame(hanam[, c("Name", "X", "Y")])[, -4]

which gives:

> head(houses)
  Name        X        Y
1 H001 20.49189 105.9274
2 H002 20.49196 105.9286
3 H003 20.49663 105.9314
4 H004 20.49191 105.9274
5 H005 20.49529 105.9303
6 H006 20.49530 105.9329

Plotting the houses

We need the OpenStreetMap package:

> if (! "sf" %in% rownames(installed.packages())) install.packages("OpenStreetMap")
> library(OpenStreetMap)
Registered S3 methods overwritten by 'ggplot2':
  method         from 
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang

Below is the code you’d need to download the tiles from 3 different types of maps:

> upperleft  <- c(20.525681, 105.905383)
> lowerright <- c(20.463343, 105.969019)
> bing <- openmap(upperleft, lowerright, type = "bing", minNumTiles = 20)
> osm  <- openmap(upperleft, lowerright, type = "osm",  minNumTiles = 20)
> esri <- openmap(upperleft, lowerright, type = "esri", minNumTiles = 20)

But the downloading would be twice as fast if you downlaod the objects directly from here:

> download.file("https://www.dropbox.com/s/26y2pgouodj4rpe/bing.rds?raw=1", "bing_file.rds")
> download.file("https://www.dropbox.com/s/rd34k3ixwthzg9g/esri.rds?raw=1", "esri_file.rds")
> download.file("https://www.dropbox.com/s/827ze3xr0orbab0/osm.rds?raw=1" ,  "osm_file.rds")
> bing <- readRDS("BING.rds")
> esri <- readRDS("ESRI.rds")
> osm  <- readRDS("OSM.rds")
> for(file in paste0(c("bing", "esri", "osm"), "_file.rds")) file.remove(file)

The function that makes the map:

> plot_hh <- function(map, points) {
+   require(sf)
+   plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)))
+   plot(map, add = TRUE)
+   plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)), add = TRUE, col = "red")
+ }

OSM map:

> plot_hh(osm, hanam)

ESRI map:

> plot_hh(esri, hanam)

BING map:

> plot_hh(bing, hanam)